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ABSTRACT 

Self-consistent solutions for triaxial mass models are highly non-unique. In general, some of these 
solutions might be dynamically unstable, making them inappropriate as descriptions of steady-state 
galaxies. Here we demonstrate for the first time the existence in triaxial galaxy models of an instability 
similar to the radial-orbit instability of spherical models. The instability manifests itself when the 
number of box orbits, with predominantly radially motions, is sufficiently large. N-body simulations 
verify that the evolution is due neither to chaotic orbits nor to departures of the model from self- 
consistency, but rather to a collective mode. The instability transforms the triaxial model into a 
more prolate, but still triaxial, configuration. Stable triaxial models are obtained when the mass 
contribution of radial orbits is reduced. The implications of our results for the shapes of dark-matter 
halos are discussed. 

Subject headings: galaxies: elliptical and lenticular, cD - stellar dynamics - methods:numerical, (cos- 
mology:) dark matter 
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1. INTRODUCTION 

Even in a collisionlcss stellar system, it is possible for 
density perturbations to grow, by inducing motions that 
reinforce the original overdensity. Such collective insta- 
bilities typically require the unperturbed motion to be 
highly correlated, and they have been most thoroughly 
studied in thin disks, which are subject to a variety of 
instabilities when sufficiently "cold." In elliptical galax- 
ies, where the stellar motions are more nearly random 
in direction, density perturbations might be expected to 
rapidly attenuate as the stars move along their respec- 
tive orbits. However it turns out that the motion in 
a variety of physically reasonable models of "hot" stel- 
lar systems is sufficient ly correlated to induce growing 
modes (as reviewed bv iMerrittl Il999f l. The two classes 
of instability that have most thoroughly been studied in 
this context are bending instabilities, which are driven 
by the centrifu gal force of stars moving across a bend 
(jToomre l fT966f): and the radial-orbit instability (hence- 
forth ROI) , which is caused by the tendency of eccentric 
orbit s to clump around a bar-like distortion (|Antonovl 
I1973L lLvnden-Belllll979h . Bending instabilities may be 
responsible for the lac k of elliptical galaxies more elon- 
gated than ~ 1 : 3 (fPolvachenko fc Shukhmanl 119771 ; 
iMerritt fc Sellwoodlll994h : simulated dark-matter halos 
also appear never to exceed this degree of elongation, pre- 
sumably because m ore flattened hal os "puff up" due to 
the instability (e.g. iBett etai] I2007D . The ROI, on the 
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other hand, can be present even in precisely spherical 
galaxies if they contain an abundance of stars on elon- 
gated orbits. It manifests itself most naturally in collapse 
simulations, which produce prolate/ triaxial bars if the 
initia l conditions are sufficiently cold (|Merritt fc Aguilarl 
1985). The ROI has also been invoked as a factor that 
regula tes the density profil es of simulated dark-matter 
halos (jBellovarv et al.ll2008f) . 

Construction of stationary models of hot stellar sys- 
tems can be difficult, and this is one reason why most of 
the studies cited above have adopted highly symmetric 
models: typically spherical, in the case of the ROI, and 
axisymmetric in the case of the bending mode studies. 
But there is no reason why such instabilities should be 
limited to spherical or axisymmetric models. Here, we 
report on a dynamical instability in triaxial models that 
closely mimics in its behavior the ROI of spherical mod- 
els. As in the spherical case, the instability leads to a 
final configuration that is close to prolate. Our simula- 
tions provide the first concrete evidence that dynamical 
instabilities may limit the permitted range of shapes of 
triaxial stellar systems, a result that may have implica- 
tions for our understanding of elliptical galaxy and dark 
halo dynamics. 

The self-consistent triaxial models on which our 
work is based were describe d in an earlier paper 
(jCapuzzo-Dolcetta et al.l l2007| ). We briefly describe 
these models in §2. The discretized models are described 
in §3, and the results of A^-body integrations in §4 and 
§5. §6 explores the dependence of the stability properties 
of the models on their orbital composition. §7 discusses 
the implications for the dynamics of elliptical galaxies 



2 



Antonini et al. 



and dark matter halos. §8 sums up. 

2. THE SELF-CONSISTENT TRIAXIAL MODELS 

The instability was discovered while testing, by TV- 
body simulations, the equilibrium characteristics of the 
triaxial galaxy models constructed in CLMV07. In this 
section we summarize the way in which the self-consistent 
orbital solutions were obtained and in the next section 
we discuss the discretized models used in the TV-body 
simulations. 

CLMV07 constructed three different self-consistent so- 
lutions of triaxial, cuspy elliptical galaxies embedded in 
triaxial dark halos. The systems differ in terms of the 
shape of the dark matter halo: (i) one halo has the same 
axis ratios as the luminous matter (1:0.86:0.7); (ii) the 
second halo has a more prolate shape (1:0.66:0.5); (hi) 
the third halo has a more oblate shape (1:0.93:0.7). Our 
choice was to study the dynamical features of the most in- 
teresting case of maximal triaxiality (i.e., the model with 
triaxiality parameter T = (a 2 — b 2 )/(a 2 — c 2 ) = 1/2), and 
to study the time evolution of the two self-consistent so- 
lutions (MODI and MODl-bis) obtained in CLMV07 for 
this case. The main difference between these two solu- 
tions was the maximum time adopted for the orbital in- 
tegrations, which was, in MODl-bis, longer (~ 5 Hubble 
times) than in MODI (~ 2 Hubble times). The models 
were constructed by means of the orbital s uperposition 
method introduced by ISchwarzsc hild (1979) which relies 
on an optimization technique. The optimization prob- 
lem consisted in minimizing the discrepancy between the 
model cell masses obtained by integration of the given 
analytical density law p(x, y, z) and those given by a lin- 
ear combination of the orbits computed in the potential 
generated by p. In our case, the two quantities to be 
independently minimized were 
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where -Bfc ;:(; im(dm) is the fraction of time that the kth or- 
bit spends in the jth cell of the luminous-matter grid 



(dark-matter grid); AL. 
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is the mass which the 



model places in the jth cell of the luminous-matter grid 
(dark-matter grid). Ck-im and Ck-dm represent the to- 
tal mass, respectively, of the luminous and dark matter 
component assigned to the fcth orbit (1 < k < n or i>). 
The basic constraints were Ckdm > and Ck-dm > 0, 
i.e., non-negative orbital weights. The departure from 
self-consistency was measured in CLMV07 by 
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where M is the average mass contained in the grid cells 
and x 2 are the quantities defined above; thus 8 represents 
the fractional rms deviation in the cell masses. 

The mass model considered in CLMV07 for the lumi- 
nous component was a triaxial generalization of Dehnen's 
(1993) spherical model with a "weak" inner cusp, p ~ 



The luminous mass density was 
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m~ = — + jj + — , < ci < bi < ai (5) 
a i b i c i 

and Mi the total luminous mass. For the dark compo- 
nent the adopted mass density was 

Pdm,0 



Pdm{m') = 



(1 + m')(l +m' 2 ) 



(G) 



/2 x- y" z" 
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u dm °dm a dm 

and Pdm,o the central dark matter density (|Burkertl 
1995). Therefore the dark component has a low-density 
core. 

In the present work we adopt the same units used in 
CLMV07: G = a; = Mi = 1. Consequently, the time 
unit is: 

[T] = G-^ai^Mr 1 / 2 (8) 

y iw^J \v=r-> ■ (9) 
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The, derived, velocity and energy units are V u = 
\J GMi I ai and E u — (GM 2 /a{), respectively. In this 
units the half mass crossing time of the system is: 
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where is the radius containing half of the model mass, 
considering the dark matter halo truncated at r = 80; 
Md is the total mass of the dark matter component. In 
the following, this time will be considered the reference 
time scale. 

3. DISCRETIZED MODELS AND THEIR PROPERTIES 

In this section we explain the methods we used to dis- 
cretize the self-consistent models described above and 
how we computed their properties for the TV-body sim- 
ulations. We also present some kinematical features of 
the models that are relevant to their stability properties. 

3.1. Discretization 

The initial conditions for the ./V-body integrations were 
set by populating the generic fcth orbit with a number 
of particles proportional to Ck and randomly choosing 
positions and velocities from the recorded data of the 
orbital integrations. In more detail: 

1 The values Ck for both the dark matter and the 
luminous matter are read from model data and the 
following quantities are evaluated: 

Nk;lm—Ck;lm/'mi m , (11) 
Nk^rn—Ck-dm/mdm , (12) 

N k =N k . M +N k;dm , (13) 

where N k is the total number of particles that pop- 
ulate the fcth orbit while N k -i m , Nk-dm are the num- 
ber of stars and dark matter particles, respectively; 
m dm and mi m are free parameters that specify the 
mass of individual star and dark matter particles. 

2 Nk } dm particles with masses mdm, and N k ,i m par- 
ticles with masses m; m , are selected with positions 
and velocities drawn uniformly and randomly (with 
respect to time) from the stored positions and ve- 
locities of the orbit integration. 
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Fig. 1. — Density profiles of the discretized model MODl-bis 
plotted versus elliptical radius m; the curves are the analytic input 
profiles. 

3 If Nk is greater than the number of available data 
in "initial data" , further positions and velocities 
are assigned using a cubic spline interpolant. 

These steps are repeated for all orbits of the self- 
consistent model. In this way the mass distributions of 
the models are adequately reproduced as long as the var- 
ious orbits are populated with a sufficiently large number 
of particles. 

The total number of objects depends on m^m and mi rn 

as 
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We chose to assign the same mass to each particle of 
the same type (luminous or dark matter), hence different 
numbers of particles are spread on orbits with different 
values of Ck- However, different values were chosen for 
the mass associated with dark and luminous components. 
We took mi — 5 x 10~ 5 and mdm = 7 x 10~ 5 which gave, 
for instance in the case of MODl-bis, a total of 166194 
particles, of which 19709 were "stars" and 146485 "dark 
matter"; the mean number of particles per orbit was 36. 
With this resolution the theoretical mass profiles were 
well reproduced as shown in Figures [T] and [2] for model 
MODl-bis. 

3.2. Determination of shapes 

The evolution of the shape of the iV-body systems was 
studied by computing the axis ratios at different dis- 
tances from the center, and also by constructing isoden- 
sity contours. 

The symmetry axes were determined from the inertia 
tensor as 



C = VTn/T„ 



( 16 ) 

where Ta are the principal moments of the inertia tensor 
and T max — max{Tn, T22, J33}. Referring to a coordi- 
nate system in which the inertia tensor is diagonal, the 
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The axis ratios of the models were compu t ed th rough the 
standard procedure descri bed by lKat3 (|1991|) ( s imilar 
met hods are describ e d by [Dubinski fc Carlbergj (|1991l ) 
and iPoon fc MerrittJ (j2004T )T To evaluate the system 
shape within a sphere of radius d, the following itera- 
tive method was used: 

1. The inertia tensor defined by particles within a 
sphere of radius d is calculated. 



2. The axis ratios are determined from equation (|16[), 

3. New axis ratios are computed considering only par- 
ticles enclosed in the ellipsoidal volume having the 
axis ratios determined in step 2. Therefore, a parti- 
cle i is included in the summations if qi < c?, where 



These three steps are iterated until the axis ratios con- 
verge. Finally, we defined a > b > c assuming c/a — 
min{C,7y, 9} and b/a the intermediate value between 
(C,rj, 8). Evaluation of the axis ratios of the discretized 
models verified the accuracy of the technique: for the so- 
lution MODl-bis we found b/a = 0.86, c/a = 0.69 for the 
luminous matter at r = 12 and b/a — 0.86, c/a = 0.72 
for the dark matter at r = 25, compared with the given 
values b/a — 0.86, c/a = 0.70 of the analytical density 
law. 

3.3. Specification of streaming motions 

X (Z)-tube orbits are here defined as those orbits hav- 
ing a non- vanishing x{z) component of the time-averaged 
angular momentum; all other, non-tube orbits (either 
box or chaotic) are defined as semi-radial orbits. Fig- 
ure [3] shows the cumulative energy distributions of the 
various orbital families in the discretized models. There 
are significant contributions from both tube and semi- 
radial orbits in both the luminous and dark components. 

A choice must be made concerning the se nse of rotation 
of pa r ticles placed in itially on tube orbits (jSchwarzschildl 
1979; Mcrritt 1980). The time-averaged density of an 
orbit is invariant to a change in sign of the initial veloc- 
ity; maximum rotation (i.e. streaming) is obtained if all 
particles on each tube orbit have the same sense of rota- 
tion, while zero mean motion is achieved by populating 
the tube orbits equally in both directions. To investi- 
gate the effects of non-zero streaming, we constructed 
two discretizations of each self-consistent solution hav- 
ing the two extreme cases of maximum and minimum net 
streaming. In Table [1] some parameters of four iV-body 
systems, sampling MODI and MODl-bis, are given. 

3.4. Model kinematics 

The kinematical features of the discretized models were 
analyzed by computing the first and the second moments 
of the stellar distribution function on a spatial grid (e.g. 
lMerrittlH98l . 
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Features of N-body models. 
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Fig. 3. — Cumulative (by mass) energy distributions of the 
various orbital families for the dark matter (DM) and the luminous 
component (LM), in the self-consistent solutions MODI (top) and 
MODl-bis (bottom). The symbols "R" , "X", "Z" and "T" denote 
the mass contributed by semi-radial, X-tube, Z-tube and tube 
orbits, respectively. 



We used different, Cartesian grids for the two different 
matter components. In the case of the dark matter grid 
the cells were cubes with sides of length 6 and all cells had 
the same size. Because of the high density concentration 
of the luminous component, we used grid cells with a 
range of sizes for the luminous matter grid. This grid 
consisted of a set of cubic cells with sizes ranging from 
0.5 near the center to 2 at greater distances. 
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Note. — The acronyms HL and LL stand for "high L" and 
"low L" where L is the absolute value of the angular momentum. 

A total of nine quantities were averaged for both mat- 
ter components in each cell: the mean velocity (Vi) and 
the six independent components of the tensor (ViVj). In 
this way the velocity dispersion tensor, 

4 = W) - <Wi>' 

could be evaluated in each cell. Then, cr^ was diago- 
nalized obtaining the three "principal" dispersions, and 
the three direction cosines giving the directions of the 
eigenvectors. 

Figure 2] shows the velocity anisotropy in the x — y 
plane for MODl-bis in the case of high angular momen- 
tum. The length of each cross arm is proportional to 
the principal value of of-. Some important features are: 
(1) a high degree of anisotropy in both components at 
all radii; (2) a nearly constant (as a function of radius) 
radial velocity dispersion of the luminous matter. On 
the other side, the radial velocity dispersion of the dark 
matter decreases strongly with radius. 

We evaluated the "anisotropy parameter" 2T r /T t 
where T r = (v*/2) and T t = (uf/2) (in full isotropy, 
2T r /T t — 1). We stress that the interpretation of this 
parameter, which is straightforward in spherical geome- 
try, is more complicated in the triaxial case. Neverthe- 
less, they give some indication of the average degree of 
velocity anisotropy. All of our discretized models yielded 
about the same values for the anisotropy parameters: 
( 2T r/Tt) dm ^ ~ 2 and {2T r /T t ) lm ~ 1.4. The high de- 
gree of "anisotropy" in the dark component - too large 
to be accounted for simply in terms of the triaxial ge- 
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ometry - suggests a strong bias toward radial motions in 
the dark matter halo, as indeed can be seen in Figure [H 

4. FIXED-POTENTIAL INTEGRATIONS 

Following l Smith fc Miller! (|1982f ) , a useful technique for 
checking whether a discrete, iV-body representation of an 
equilibrium model was correctly constructed is to inte- 
grate the discretized initial conditions in the fixed poten- 
tial of the analytic model and observe whether there is 
any change in the spatial distribution of the N particles. 
This procedure is also a powerful way to constrain the na- 
ture of any evolution that is observed in the full iV-body 
integrations. Actually, if the shape of a model changes 
with time, this could be due either to chaotic evolution 
of individual orbits (mixing) , or to collective modes (dy- 
namical instability). But the former mechanism will be 
active even when the potential is fixed, whereas a collec- 
tive mode requires an evolving potential. 

Accordingly, we integrated the orbits of the discretized 
models in the analytic potential corresponding to the 
smooth mass distributions (equations l2l2| . Each particle 
was advanced, up to 30 crossing times, using a 7/8 or - 
der Runge-Kutta algorithm described bv lFehlberel (|1968l ) 
with a variable time step, in order to keep the relative er- 
ror per step in energy less than a specified value (10 -8 ). 
Since each orbit is independent, this operation is easily 
parallelized. The simulation required ~ 5 CPU hours 
total for 166000 particles. Registration of the particle 
positions and velocities were made at fixed intervals of 
time, the same for all particles. The duration of 30t cross 
is longer than the time over which the instability mani- 
fests itself in the full iV-body simulations (see below). 

Figure [5] shows that no significant evolution of the axis 
ratios is observed in the fixed potential integrations. Also 
the contours of the projected density for both systems 
remain approximately unchanged until the end of the in- 
tegrations. Actually, the relative variations of the axis 
ratios with respect to their initial values are within 4% 
for the luminous component and 1% for the dark mat- 
ter. Given the unavoidable noise in the computation of 
the axis ratios, such variations are irrelevant; the larger 
variations in the luminous matter are probably a con- 
sequence of the higher noise due the lower number of 
particles. These results allow us to conclude that the 
initial conditions were correctly generated, and also that 
any strong global shape deformations in the full ./V-body 
simulations must be a manifestation 

of a dynamical instability, and not chaotic mixing of 
individual orbits. 

5. Af-BODY INTEGRATIONS 

5.1. N -Body code 

The full iV-body integ rations were carried out us- 
ing th e TreeATD code of iMiocchi fc Capuzzo-Dolcettal 
(2002), which is a parallel code that uses a tree algorithm 
for the gravitational force evaluation, and an individual 
time stepping for the leap-frog integrator. TreeATD needs 
three input parameters that influence the speed and ac- 
curacy of a simulation: the opening angle 9, the smooth- 
ing length e, and the maximum allowed time step At. 
We set e = 0.05 and e = 0.1 in the case of model LL},i S , 
At, = 0.07 and 9 = 0.7. These values were chosen in or- 
der to conserve energy within 0.05% over the full course 
of the integrations. Simulations were performed using 



8 nodes of gravitySimulator, a 32-node cluster at the 
Rochester Institute of Technology. 



5.2. The instability 

The ./V-body integrations revealed that both MODI 
and MODl-bis represent unstable equilibria, in the sense 
that their axis ratios evolve significantly. Lagrangian 
radii for both models showed essentially no evolution, 
indicating that the instability affects only the shapes of 
the models and not the global concentration of matter. 

Figure [6] illustrates the change of the axis ratios up to 
t = 711.2 (40 crossing times) for the maximum-streaming 
models HL and HL^is- Strong deformations appear ev- 
idently in both the luminous and dark components after 
just two crossing times. The initial evolution is toward 
a more spherical shape; the duration of this first phase 
is different in the two components, being longer for the 
dark component than for the luminous component. Fi- 
nal shapes are nearly prolate, with 0.7 < b/a < 0.77 and 
0.65 < c/a < 0.7 for the dark halo and 0.69 < b/a < 0.71 
and 0.59 < c/a < 0.62 for the luminous matter. The fi- 
nal contours of the projected density for model HL are 
shown in Figure [7] After ~ 25 crossing times (dark mat- 
ter) and 18 crossing times (luminous matter), the insta- 
bility appears to have run its course, but the system still 
exhibits a slow figure rotation, as shown in Figure [5] 

In order to understand better the influence of model 
rotation on the dynamical evolution, we performed a sec- 
ond set of simulations for the systems LL and LLu s , 
which correspond to the case of minimum angular mo- 
mentum. The right columns of Figure [6] display the evo- 
lution of the axis ratios for models LL and LLu s , while 
Figure O shows the contours of the projected density for 
both the mass components at the final time for model LL. 
The results confirm that figure rotation is suppressed in 
these cases; nevertheless, dynamical evolution still leads 
to the formation of a sort of bar. A comparison between 
the evolution in the cases of high and low L suggests 
that a dynamical instability occurs at the beginning and 
manifests itself completely in ~ 20 crossing times. Af- 
ter this time, figure rotation is still present for HL and 
HLbi S while the LL and LLu s systems conserve both 
their shape and their orientation in space. 

Based on Figure [6l the change in shape occurs sooner 
in the case of low angular momentum; for models LL 
and LLbi S , the "stable" phase begins at ~ 8t cross in the 
luminous component and at ~ 18t cross in the dark com- 
ponent. In addition, in the absence of rotation, the final 
elongation in the dark matter is greater: 0.64 < b/a < 
0.69 and 0.6 < c/a < 0.67. Rotational motion has the 
effect of breaking slightly the axisymmetry reached by 
the model. 

This evolution is strongly reminiscent of the well- 
kno wn ROI seen i n radially-anisotropic, spherical mod- 
els (jMerrittl I1999T ). In spherical models, the instabil- 
ity causes a bar to form with some random orientation, 
determined by the precise spectrum of density inhomo- 
geneities in the initial model. In the triaxial case, the 
initial conditions arc already bar-like, but the instabil- 
ity chooses a new bar-like distortion to grow. As in the 
spherical case, the final configuration after the instability 
has run its course is close to prolate. 
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Fig. 4. — Velocity dispersions of MODI — bis in the case of luminous matter (left) and for the dark matter (right). The components of the 
velocity dispersions were stored in each cell of a grid divided as described in the text. The length of the arms of each cross is proportional 
to the corresponding principal dispersion. 
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Fig. 5. — Upper panels: evolution of the axis ratios for the so- 
lution MODl-bis in the fixed-potential integrations. Lower panels: 
evolution of the axis ratios for the solution MODI. The axis ratios 
are evaluated at r = 8 for the luminous matter and at r = 25 for 
the dark component. The times are scaled to the internal crossing 
time. 



6. 



DEPENDENCE OF THE INSTABILITY ON THE 
ORBITAL COMPOSITION 



In spherical models, the ROI is associated with a pre- 
dominance of eccentric orbits. The instability growth 
rate can be reduced to zero by changing the orbital com- 
position toward more isotropic or tangentially-biased so- 
lutions; this is always possible since there are many dis- 
tribution functions f(E, J) that correspond to the same 
density profile p(r). Likewise, in the triaxial geometry, 



the Schwarzschild method can yield a variety of orbital 
solutions consistent with a specified mass model, and 
these different solutions will generally have different sta- 
bility properties. Here we show that the instability de- 
scribed above in the triaxial models can indeed be effec- 
tively suppressed by reducing the number of semi-radial 
(box) orbits in the self-consistent solutions. This result 
reinforces our hypothesis that the instability is intrin- 
sically similar in character to the ROI, and is also of 
physical significance for real galaxies, as discussed in §7. 

6.1. Minimizing the contribution from semi-radial orbits 

To investigate the hypothesis of ROI, new orbital solu- 
tions constraining the number of semi-radial orbits were 
constructed, discretized and evolved forward in time as 
iV-body sys tems. 

Following iPoon fc Merrittl (|2004f ) , the relative contri- 
butions of different orbits to the self-consistent solutions 
was varied by adding a penalty function to equations JT]) 
and @, which became 




k;lm^k^j\hn 



(19) 



and 



Xdn 



j;dm 




E 

k=l 



\drn ■dm 



(20) 



respectively. Here, 



k\lm{dm) 



is a penalty associated 



with the fcth orbit of the luminous (dark) component; 
as Wfc ; i m (dm) increases, the mass contribution C 'k-im{dm) 
of the fcth orbit in the model decreases. (We remark that 
the role of our penalty function is that of an ad hoc nu- 
merical device and does not have any particular physical 
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Fig. 8. — Evolution of the isodcnsity contours in the y — z plane for the model HLi, is ; the two horizontal set of panels on the top display 
the configuration at 0, 4, 8, 12, 16t cross , then, the two bottom set of panels refer to 20, 24, 28, 32, 36t cross At each time, the upper panels 
represent the luminous component while the dark matter is shown into the lower panels. In the case of the luminous matter the linear size 
of each box is 22 while for the dark matter it is 85. The arrows represent the eigenvectors of the inertia tensor with the highest value of 
the projection on the y — z plane (in most of the plots, the other eigenvector is approximately along the line of sight). The rotation can 
clearly be seen after 20 crossing times; by this time the instability appears to have run its course. 



t 1 1 1 — I I 1 1 1 1 r 




Fig. 9. — Contours of the projected density at 40 crossing times for model LL. Left panels: luminous component; right panels: dark 
matter component. 
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Fig. 10. — Cumulative energy distributions (by mass) of the 
various orbital families for different values of Wr. The symbols 
R, X, Z and T denote semi-radial, Jf-tube, Z-tube and all tube 
orbits, respectively. 

meaning.) We chose Wk dm = Wk-dm = for the tube 
orbits and Wkdm = W R . dm > and W k - dm = W R - dm > 
for the semi-radial orbits. The optimization problem 
represented by equations (|19[) and (|2T)]) was solved us- 
ing the NAG routine E04NCF, which implements an ef- 
ficient method to solve solves linearly constrained linear 
least-squares pr oblems and convex qu adratic program- 
ming problems (|Stoerl Il97it lGill|[l9M ). The new solu- 
tions were found using the full orbital library correspond- 
ing to model MODl-bis. 

As shown in Figure fTUl setting both W R; d m > W R . dm > 
does in fact increase the number of tube orbits in the 
solutions at the expense of the semi-radial orbits. The 
error in the cell masses increases at increasing W R , as 
shown in Figure 111! How ever, it is well known that (e.g. 
iMerritt fc Frid man 199(3) 5 the value of 5 alone is not able 
to judge the degree of self-consistency of an orbital solu- 
tion. Therefore, our new solutions might still represent 
reasonable equilibria even if the quality of the fit to the 
cell masses is worse than that of the solutions found in 
CLMV07. We return to this point below. 

6.2. Discretized models and their kinematical properties 

We examined the properties of discretized models for 
four choices of the penalty parameters: 

(i) W R - dm = W Rdm = 50 (model Nl); 

(ii) W R , dm = W R , lm = 5 (model N2); 

(iii W R - dm = 5 x 1CT 3 and W Rdm = 5 x 10~ 6 (model 
iV3); 

(iv) W R , dm = 50 and W Rdm = 5 x lO" 6 (model N4). 
./V-body realizations were generated as explained in §3; 
the sense of circulation of the tube orbits was chosen 
randomly. Each model used 2 x 10 4 luminous particles 
and 1.5 x 10 5 dark matter particles. Table [2] gives values 
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Fig. 11. — Departure from self-consistency (<5) as a function of 
the "penalty" (Wr). 



TABLE 2 

Anisotropy parameters of the new models. 



MODEL 


w mm 


w R . dm 


[2T r /T t \i m 


[2T r /T t \ dm 


Nl 


50 


50 


0.512 


1.175 


N2 


5 


5 


0.784 


1.335 


N3 


5 x 10~ 6 


5 x 10~ 3 


1.220 


1.754 


N4 


5 x 10~ 6 


50 


1.230 


1.174 



of the anisotropy parameters. As expected, (2T r /T t )i m 
and (2T r /T t ) dm are decreasing functions of W R . 

Figure fT2l displays the velocity dispersions of these sys- 
tems in the x — y plane, clarifying the relation between 
W R and the velocity anisotropy: when W Rdm =5 or 50, 
the tangential velocity dispersion is higher than the ra- 
dial dispersion in the luminous matter. For W Rdm = 
5 x 10^ 6 , we found a r w er t . In the case of the dark 
component, when W Rtdm = 5 or 50 the halo is nearly 
isotropic; in these systems the principal axes of the ve- 
locity ellipsoid tend to lose their radial alignment. In 
the case W R - dm — 5 x 10~ 3 , the dark matter becomes 
strongly anisotropic. 



6.3. N-body simulations 

The new ./V-body systems were evolved for 20 crossing 
times, approximately the time scale of the instability to 
grow, as seen before.. Figure [T3l shows the evolution of 
the model axis lengths. Model A^3 is clearly dynamically 
unstable, evolving into a prolate configuration. By con- 
trast, for models Nl, N2 and NA, the instability was 
either absent or much suppressed. Figures [14] and [15] 
show contours of the projected density at t = 20t cross 
for models Nl and iV4 respectively; in these cases, the 
final configuration looks very similar to the initial one. 

The correlation between bar formation and the value 
of the radial velocity dispersion is a clear sign that the 
dynamical instability discussed in §4.2 can be identified 
with the ROI. Furthermore, the behavior of model NA 
suggests that the instability disappears when the dark 
halo is made isotropic, i.e., the instability derives mainly 
from anisotropy in the dark component. 
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IV, 



We can draw the following conclusions from these in- 
tegrations: 

1. The dynamical instability disappears when the frac- 
tion of semi-radial orbits in the models is decreased. 

2. The stability properties of these models are deter- 
mined primarily by the kinematics of the dark matter 
halo. In particular, stable configurations are obtained 
when {2T r /T t ) dm < 1.4. 

7. DISCUSSION 

Large-scale simulations of structure formation have 
shown that dark matter halos, during their evolu- 



tion, develop universal properties, such as a char- 
acteristic density profile, a power-law dependence 
of phase space density on radius, a linear rela- 
tion between the velocity anisotropy and the den- 
sity slope {(3 — 7 relation) and a particular distribu- 
tion of shapes and spins (Dubinski & Carlberg 119911; 
Warren et all Il992t [Navarro et alj 119961: iMoore et al l 



1999t iTavlor fc Navarrol 120011: iHansen fc Moord 12006. 
Vario us authors (e.g. ISver fc White! 119981 : iDekel et al.1 
2003) have argued that dynamical processes during merg- 
ers may be responsible for the apparent universality of 
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Contours of the projected 



density at 20 crossing times in model iVl for both the luminous matter ) and the dark matter 
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Fig. 16. — (a) Axis ratios of the dark matter halos of our evolved unstable models. LL indicates models with low angular momentum 
while HL refers to rotating models. IM refers to the axial ratios of the initial model, (b) Mean axis ratios of dark matter halos in simulations 
of structure formation in ACDM and AWDM cosmologies, (c) Shapes of eq uilibrium JV-body mode ls formed by isolated cold collapse. 
Models that were formed from rotating initial conditions are labelled ROT. The is fc Spurzcm (1999) gi ve, fo r a set of 9 dissipationless 
collapse simulations, the axis ratios of the 50% and 10% of the most bound particles. For lBeUovarv et al.l (120081 ) we show two values of the 
final axis ratios of their initially coldest model: one refers to the innermost regions (core) and the other one to a larger distance, roughly 
the virial radius. 
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thes e relations. On the o ther hand, iHuss et all (|1999f ) 
and lWang fc White! (|2008| ) have pointed out that many 
of these properties are also reproduced in a universe 
where halos form via monolithic collapse, suggesting that 
mergers are not essential for establishing the universal 
relations. 

Some of the regularity in dark matter halo proper- 
ties may be due to dynamical instabilities, which limit 
the range of allowed equilibrium states irrespective of 
how the halos formed. A well-established example is 
the effect of bending instabilities on the shapes of hot 
stellar systems: major to minor axis ratios are lim- 
ited to ~ 3 : 1 for both oblate and prolate systems 
dMerritt fc Sellwoodl fl99l . This is a plausible expla- 
nation for th e lack of elliptical galaxies flatt er than Hub- 
ble type E7 ijFridman fc Polvachenkolll984l ), and is also 
consistent with the maxim um elongations found for sim- 
ulated dark matter h alos (|Bullocki 120021 : lAllgood et al.1 
120061: IBett eHiIll2l)07l) . 

The role of the ROI in establishing such "universal" 
characteristics is less clear. The ROI arises naturally 
in halos formed via monolithic collapse, causing oth- 
erwise spherical systems to settle into prolate/triaxial 
shapes dMerritt fc Aguilarll985t[Aguilar fc Merritll99d ; 
IHuss et al.lll999f) . The instability also reduces the depen- 
dence of the final concentratio n on the initial "temper- 
ature" of the collapsing cloud dMerritt fc Aguilari Il985t 
iBarnes et al.ll2005t iBellovarv et al.ll2008l ). But formation 
via mergers is qualitatively different from collapse; and 
since a (spherical) model can always be rendered stable 
by making its velocity distribution sufficiently isotropic, 
the role that the ROI plays in determining the structure 
of dark matter halos is likely to depend somewhat on the 
details of the halo formation process. 

We nevertheless note that the halos formed in hier- 
archical cosmologies ten d to exhibit radia l ly-anisotropic 
envelopes, a r /&t ~ 1-5 dColiri et all 120001 : IWojtak et al.1 
l2005t iNavarro et al.l l2008h . and that these anisotropies 
are similar to those of unstable models formed via col- 
lapse, both before and after the ROI has run its course 
dBellovarv et al.l 120081 ) , and consistent with the values 
that render our two-component models unstable. So it 
is plausible that the ROI or something similar is active 
during the hierarchical formation of halos. 

Figure [TBI presents a weak test of this idea. Axis ratios 
of our unstable halo models (panel o) are compared with 
those of dark matter halos formed in various cosmological 
simulations (panel b) and with TV-body models formed 
via simulations of isolated collapse (panel c). As noted 
above, the instability has the effect of making our ini- 
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tially triaxial models more prolate and more elongated. 
In the absence of rotation (LL), the final shapes are more 
prolate than for typical cosmological halos. However in 
the models with streaming motions (HL), the final axis 
ratios are essentially identical to the average values found 
in the cosmological simulations. A similar conclusion can 
be drawn from the isolated collapses in panel (c), which 
also tend to be more triaxial (i.e. less prolate) when 
rotation is present. 

We stress again that our unstable models could be 
rendered stable by selecting different sets of orbits, in 
the same way that unstable spherical models can be 
stabilized by making their velocity distributions more 
isotropic. The role of the ROI in structuring dark mat- 
ter halos must therefore depend somewhat on the orbital 
composition of halos formed in the cosmological simula- 
tions. 

8. CONCLUSIONS 

We explored the stability properties of two, self- 
consistent models of triaxial galaxies embedded in tri- 
axial dark matter halos. Our results can be summarized 
as follows. 

1. Both models were found to be dynamically unsta- 
ble, evolving toward more prolate shapes on a time scale 
of ~ 20 crossing times. Final shapes were approximately 
prolate in both components, with short-to-long axis ra- 
tios of - 0.6 - 0.7. 

2. The evolution was shown not to be due to errors in 
construction of the equilibrium models, nor to diffusion 
of chaotic orbits, but rather to a collective mode. On 
this basis we identified the instability with the ROI of 
spherical models. 

3. Including streaming motions in the initial models 
leads to final configurations that are more triaxial than 
when rotation is absent. These final shapes are very sim- 
ilar to the mean shapes of dark matter halos formed in 
hierarchical merger simulations. 

4. When the number of box-like orbits is reduced 
below a certain threshold the dynamical instability 
disappears. The presence or absence of the instability 
is most strongly affected by the number of box-like 
orbits in the dark matter halo; stable configurations are 
obtained when (2T r /T t ) dm < 1.4. 
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